On-orbit servicing is a vital method to extend the lifetime and enhance the performance of spacecraft.^[1] One of the challenging techniques is the autonomous proximity. For cooperative spacecraft, several space missions^[2,3] have demonstrated the technique in orbit. However, for a non-cooperative object such as a tumbling spacecraft, it is still an undemonstrated risky operation.^[4] With the rotation of the target, the chaser needs to approach to and align with the target to prepare for the final docking, capture, or photographing. This involves the precise position and attitude synchronization simultaneously, and the accurate relative motion modeling as well as nonlinear control design for the highly nonlinear kinematics and dynamics is crucially important.^[5,6]

Traditionally, rigid-body dynamics can be represented by translation and rotation about the center of mass (CM).^[7–11] However, during the final phase of proximity, the relative position and pointing of the spacecraft’s sensors or devices are usually worth paying more attention to, as they need to satisfy some specified constraints, such as the field-of-view constraint,^[12] collision avoidance constraint, and plume impingement constraint.^[13] So, it is advisable to model and control the motion of these sensors and devices directly. As these instruments are usually equipped on the surface of the spacecraft, not coincident with the CM, their motion can be affected by the translation of the CM and the rotation about the CM. Thus, the translation and rotation are kinematically coupled. This coupling was first addressed by Segal and Gurfil^[14] under the attitude synchronization assumption. Six years later, Lee et al.^[15] abandoned the assumption and derived the kinematically coupled nonlinear relative translational and rotational motion model in matrix form. As the translation and rotation dynamics were set up in the local-vertical local-horizon (LVLH) frame and chaser’s body-fixed frame, respectively, the kinematical coupling led the whole model to be quite complex.

As the refueling of spacecraft is costly and risky, the fuel consumption for the controllers is a vital performance index, besides the control accuracy. While the design of optimal controllers for nonlinear systems is a very challenging topic, one of the outstanding techniques is based on the state-dependent Riccati equation (SDRE) method,^[16] which transforms the nonlinear system into a linear-like structure with state-dependent coefficient (SDC) matrices, and minimizes the quadratic index. Although many important results^[17–19] for regulator problems have been obtained, the optimal tracking control based on the SDRE method is still a hot topic. In Refs. [9], [11], and[20], the SDRE was implemented as an integral servomechanism to achieve tracking. Strano and Terzo^[21] proposed an SDRE-based tracking controller consisting of the feedback gain from SDRE and the feedforward term specified with the hydraulic actuation system. Recently, Batmani et al.^[22,23] pointed out that the shortage hindering the solution of the optimal tracking problem with the SDRE method is that the quadratic cost function is only valid for the desired trajectories generated by an asymptotically stable system. To eliminate the defect, a discounted cost function was adopted to transform the tracking problem into a regulator problem, deriving the suboptimal control law with SDRE.

Considering the aforementioned facts, in this paper we investigate the optimal approach to a tumbling target. Regarding that the devices or sensors are usually not consistent with the CM, a kinematically coupled 6 degrees-of-freedom(DOF) model is set up to describe the motion of the off-CM feature point on individual spacecraft. Different from Refs. [14] and [15], the attitude dynamics is uniformly established in the LVLH frame, simplifying the representation of the kinematical coupling terms. With the rotation of the target, the position and attitude trajectories for the chaser to track oscillate with time and never approach to the origin. To perform the optimal tracking, a feedback–feedforward controller is designed, with the feedforward part providing the necessary input for following the trajectory and the feedback term based on SDRE stabilizing the tracking error dynamics. The optimality of the proposed controller is analyzed, and a comparison between the proposed method and the method in Refs. [22] and[23] is made.

The rest of this paper is organized as follows. In Section 2, we give the problem formulation and equations of motion. The control methods are presented in Section 3, and the SDC matrices are also given in this section. The numerical simulations are provided in Section 4. Finally, some conclusions are drawn from the present study in Section 5.

To set up the 6-DOF model for the spacecraft components (sensors or devices), we assume that a virtual reference spacecraft moves on the Kepler orbit. As shown in Fig. 1, the Earth-centered inertial (ECI) frame is represented by the geocentric-equatorial frame {I} = {I_x,I_y,I_z}, and the LVLH frame, {L} = {L_x, L_L,L_L}, is built on the CM of the reference spacecraft, with L_L being a unit vector directed from the spacecraft radially outward, L_z normal to the orbit plane, and L_L completing the right-hand coordinate system. The rigid-body spacecraft flies close to the reference spacecraft, and the body-fixed frame is defined as {B} = {B_x, L_B, L_B}. A feature point (component) on the spacecraft (that does not coincide with the CM) is denoted by Pⁱ, and Pⁱ is its position vector in {B}. The ρ and ρⁱ are the position vectors of the CM and the feature point Pⁱ in the LVLH frame, respectively.

Fig. 1.

	Figure Option ViewDownloadNew Window
	Fig. 1. (color online) Rigid-body spacecraft with the reference LVLH frame and body-fixed frame.

2.1. 6-DOF relative motion model

From Fig. 1, one can note that the following geometric relationship holds:

and with the rotation of the LVLH frame {L} and body-fixed frame {B}, a relative angular velocity {ω} is generated aswhere ω_B and ω_L are their own angular velocities, respectively, and is the rotation matrix from frame {B} to frame {L}.

Taking the derivative of Eqs. (1) and (2) with respect to time, one can derive the kinematic and dynamic equations of the 6-DOF rigid-body spacecraft for the feature point as follows:kinematics:

dynamics:where q = [q₀, q₁, q₂, q₃]^T is the quaternion from frame {B} to frame {L}; and are their own angular accelerator vectors, respectively; is the inverse matrix of , and^[24]

The translational dynamics of the spacecraft CM, , and the rotational dynamics of the spacecraft, , can be specified respectively by^[24]

where x, y, and z are the components of ρ in the frame {L}, r₀ and f are the current orbit radius and true anomaly of the virtual reference spacecraft, respectively, μ is the gravitational parameter, and F = [F_x,F_y,F_z]^T is the control acceleration in the frame {L}, J and M are the moment of inertia tensor of the rigid body and the total external torque, respectively.

Substituting Eqs. (6) and (7) into Eq. (4) with Eqs. (1)–(3), one can obtain

where [x_i,y_i,z_i]^T and are the components of the vectors ρⁱ and Pⁱ in the frame {L}, respectively, , and .

2.2. Transformation of body-fixed frame

During the final phase of approaching to the tumbling target, the boresight of the chaser’s feature device needs to be kept aiming at the target’s feature device to perform docking, capturing, and photographing. Traditionally, the attitude alignment is achieved by making the two spacecraft’s body-fixed frames coincide with each other. However, the feature components may fail to aim at each other with attitude synchronization, if the body-fixed frames are not defined in the proper way. As shown in Fig. 2, the chaser needs to take photographs for the feature zone on the tumbling target; however, the camera cannot aim at the feature zone when the chaser body-fixed frame {B_c} = {B_{c_x}, B_{c_y}, B_{c_z}} is coincident with the target frame {B_t1 = {B_{t_x1}, B_{t_y1}, B_{t_z1}, which is defined in its own principal axes of inertia. To eliminate the defect, the target frame needs to be redefined in a proper way, {B_t2} = {B_{t_x2}, B_{t_y2}, B_{t_z2}}. The quaternion for the target to transform its body-fixed frame into the new frame can be specified as follows.

Fig. 2.

	Figure Option ViewDownloadNew Window
	Fig. 2. (color online) Contradiction between attitude synchronization and component alignment.

Assuming that the unit normal vector of the feature point on target is e_t in {B_t1}, and the unit normal vector of the feature point on the chaser e_c, the rotation axis and the angle for the target to transform its body-fixed frame {B_t1} into the new frame {B_t2} are

The corresponding quaternion iswhere , , and are the components of .

Then the inertia matrix of the target also has a change. Based on the definition of the inertia matrix, one can calculate the new inertial matrix in the new body-fixed frame from the following expression:

where p is the position vector of the infinitesimal body element dm in the frame {B_t1}, I_n is the n-dimension unit matrix, is the rotation matrix, transforming a vector from the frame {B_t1} to the frame {B_t2}, and J is the original inertia matrix in {B_t1}.

To perform the final approach, the chaser’s feature component needs to track the target’s with boresight alignment. Considering the fuel consumption, the index of control accuracy and fuel consumption is written in the quadratic form and the SDRE technique is utilized to design the suboptimal tracking controllers. Since the reference system is not asymptotically stable, the traditional quadratic cost function used in the SDRE technique is not valid.^[23] To eliminate this defect, we add a desired control signal into the index function, and propose a feedback–feedforward control structure. The desired control signal provides the necessary input to follow the reference system, acting as the feedforward term, and the feedback term is generated with SDRE. The idea is from Ref. [21], but a general method to derive the feedforward term is proposed. In addition, the existing method in Refs. [22] and [23] is also introduced for comparison. To make the difference, we mark the new method as the desired control method, and the method in Refs. [22] and [23] as the discounted cost method.

3.1. Desired control method

3.1.1. Controller structure

Consider the general nonlinear plant which is full-state observable and affine in the input

where x ∈ ℝ^{n × 1} and u ∈ ℝ^{p × 1} are the state and control input, respectively. The system is expected to track the motion trajectory given bywhere z ∈ ℝ^{n × 1} has the same dimension as x.

Under the assumption of f(0) = 0, the nonlinear function f(x) can be written in the SDC form

and the pair {A(x), b(x)} is assumed to be pointwise stabilized in the linear sense for all x in the domain of interest.

To make system (13) follow the desired state trajectory z, the performance index to be minimized is defined as^[21]

where Q and R are the positive-semidefinite and positive-definite symmetric matrices penalizing the deviations of state and control, respectively. The term u_d is the desired control signal meeting the following constraint:Then feedforward gain u_d can be specified asBesides the feedforward gain u_d, the feedback gain is generated with the SDRE technique, and the whole control law is given by^[21] where P(x) is the solution of the algebraic state-dependent Riccati equation

3.1.2. Optimality analysis

To analyze the optimality of the proposed controller, we first transform the optimal tracking problem into the optimal regulation problem. Denoting e = x − z, and w = u − u_d, and combining with Eqs. (15) and (17), the error dynamics for the tracking problems (13) and (14) can be set up as follows:

The performance index function (16) turns out to be

Bearing in mind that z can be specified with the reference dynamics (14) and the given initial state z₀, A(x) = A(e + z) and b(x) = b(e + z) each are a function of e. In addition, the proposed controller (19) under the error dynamics becomes

The Hamilton function for the infinite-horizon optimal regulator problem (22) and (23) is

where λ is the introduced costate state. Based on the Pontryagin maximum principle,^[25] the optimal solution for the regulator problem isand the following conditions along the optimal trajectory must be satisfied:where the upper equation is the Hamilton–Jacobi–Bellman(HJB) equation, and the lower three equations are the first-order necessary conditions for optimality.

As A(x) = A(e + z) and b(x) = b(e + z) each are treated as a function of e, and P(x) is the solution of the Riccati equation (20), the controller (24) is just the SDRE regulator for the error dynamics (22) and the index function (23). For the SDRE regulator, Mracek and Cloutier^[26,27] have addressed its optimality by considering the conditions (27). The main optimality results are presented as follows.

3.2. Discounted cost method

The discounted cost method exerts a discount factor on the performance index as follows:^[22,23]

where γ > 0 is the discounted factor. Defining the new state and control, and , the cost function turns to beand the corresponding weights matrices areWriting the systems (13) and (14) in the quasi-linear forms with f(x) = A(x)x and g(z) = G(z)z, the two systems can be augmented asThen the tracking problem is converted into a standard nonlinear optimal regulator problem. Based on SDRE, the control law iswhere is the solution of the Riccati equationand the suboptimal tracking controller for the origin tracking problem is given by^[22,23]

Further, the matrix can be parted as

where P₁₁ (x) ∈ ℝ^{n × n}, and the control law can be written aswhere P₁₁(x) and P₁₂(x,z) are respectively the solutions of the following Riccati and Lyapunov equations:

Then we can make a comparison between the two methods. Both methods transform the optimal tracking problem into a regulator problem, and can be seen as the derivatives of the SDRE method. The two controllers (19) and (36) are both composed of two terms. The first is the feedback part, and the same-dimension Riccati equations (20) and (37) have to be solved respectively at each time instant. The second is the feedforward part, but the calculations are different. For the desired control method, the feedforward part u_d can be derived directly from Eq. (18). While for the discounted cost method, the Lyapunov equation (38) has to be solved. So the desired control method is a little more efficient.

3.3. SDC parameterization

To fulfill the two methods, the 6-DOF motion model in Subsection 2.1 has to be rewritten in a linear-like structure with SDC matrices. Before being parameterized, the term in the attitude dynamics (9) needs to be converted as follows:

where the operator ( )^× converts three-dimensional vector η into the skew-symmetric matrix defined byWith , the term can be rewritten in linear-like form dependent on q Then, equation (9) is rewritten as

Denoting and taking the SDC parameterization of Eqs. (3), (8), and (42), the particular A(x) is selected as

where I_{n × n} is an n × n identity matrix, I_{m × n} is an m × n zero matrix, −10⁻⁴I_{4 × 4} is added to the quaternion kinematics for the Riccati equation to have a stable solution,^[11] and

The control distribution matrix is

For the desired control method, only the SDC parameterization A(x) and b(x) are calculated using Eqs. (43) and (51) at every moment. While for the discounted cost method, the reference system has also to be written in the quasi-linear form, and equation (43) is utilized again to obtain G(z).

In this section, simulation results for the scenario that the chaser needs to photograph a particular area on the tumbling target are given. The reference LVLH frame is set up on the CM of the tumbling object. The initial orbital elements of the target are listed in Table 1. The moment of inertia tensor of the target is J_t = diag([2000,2500,2000]) kg · m². To take photographs of the particular area, the virtual point P^t = 20 × [0.5567,0.6634,0.5]^T m in the target’s body-fixed frame is selected to set the chaser’s camera, and the unit normal vector of the particular area is parallel to P^t. The inertia matrix of the chaser is J_c = 500 × diag([1,1,1]) kg · m². The position of the camera in the chaser’s body-fixed frame is P^c = [0,0,0.5]^T m, and the camera’s boresight vector is parallel to P^c. The initial states for the scenario are listed in Table 2.

Table 1.

Table 1.

Table 1. Orbital elements. . Orbital elements Values Semi-major axis/km 7170 Eccentricity 0.05 Inclination/(°) 15 Right ascension of ascending node/(°) 0 Argument of perigee/(°) 340 True anomaly/(°) 20

Table 1. Orbital elements. .

Table 2.

Table 2.

Table 2. Initial states of target and chaser. . Initial conditions Target spacecraft Chaser spacecraft Relative position of CM/m ρt = [0, 0, 0]T ρc = [30, 30, 30]T Relative velocity of CM/m·s−1 Quaternion qt = [1,0,0,0]T qc = [0.5,0.5,0.5,0.5]T Relative angular velocity/rad·s−1 ωt = [0.018,0.018,0.016]T ωc = [0.035,0.105,0.069]T Relative position of the feature point/m Relative velocity of the feature point/m·s−1

Table 2. Initial states of target and chaser. .

It can be seen that the chaser camera can never aim at the target’s particular area when their body-fixed frames are synchronal. To eliminate the defect, an alternative body-fixed frame for the target is chosen. Based on Subsection 2.2, the quaternion from the original body-fixed frame to the new frame is . The relevant transformation matrix and the new inertia matrix are

For the given initial states of the target, figure 3 shows the position trajectory of the target’s feature point in the LVLH frame, which is also the desired trajectory for the chaser. The CM of the target keeps constant at the origin, so the feature point moves on the sphere with the center at the origin and a radius of 20 m. Figure 4 presents the attitude trajectory of the reference system. It can be seen that the quaternion’s components oscillate with time and never approach to zero.

Fig. 3.

	Figure Option ViewDownloadNew Window
	Fig. 3. (color online) Position trajectory of the target’s feature point in the LVLH frame.

Fig. 4.

	Figure Option ViewDownloadNew Window
	Fig. 4. (color online) Attitude trajectory of the target’s body-fixed frame.

To track the target’s position and attitude trajectory, the two controllers (19) and (36) are executed respectively. Figures 5 and 6 show the norms of position and velocity tracking errors respectively with the two methods. The full curve represents the results from the desired control method, and the dotted curve refers to the discounted cost method. In Fig. 5, the decreasing trends of position error for both controllers are similar, but the desired control method achieves a quicker convergence. The position errors during the last 100 s are less than 10⁻⁴ m for the desired control method and 0.2 m for the discounted cost method, respectively. Figure 6 shows that the two methods also have similar characteristics. The linear velocity errors during the last 100 s for the two methods are 2 × 10⁻⁶ m/s and 0.007 m/s, respectively. It can be seen that the desired control method achieves better tracking accuracy for the translation motion.

Fig. 5.

	Figure Option ViewDownloadNew Window
	Fig. 5. (color online) Norms of position tracking errors.

Fig. 6.

	Figure Option ViewDownloadNew Window
	Fig. 6. (color online) Norms of linear velocity tracking errors.

Figure 7 gives the roll, pitch, and yaw in Euler angles expressed in the 3–2–1 rotation sequence from the chaser body-fixed frame to the target body-fixed frame. The Euler angles of two controllers both converge quickly within the initial 200 s, and to less than 0.15° finally, and the desired control method’s accuracy is higher. Figure 8 shows the tracking errors of the relative angular velocity with the two methods. Also steady errors show that the desired control method is superior to the discounted cost method.

Fig. 7.

	Figure Option ViewDownloadNew Window
	Fig. 7. (color online) Euler angles from the chaser’s body-fixed frame to the target’s.

Fig. 8.

	Figure Option ViewDownloadNew Window
	Fig. 8. (color online) Tracking errors of relative angular velocity.

Figures 9 and 10 depict the control accelerations and torques respectively for the two controllers. By integrating the norms of the control inputs, we can derive that the total fuel consumptions in control acceleration are 9.01 m/s for the desired control method and 8.78 m/s for the discounted cost method, and the total consumptions in control torque are 184.48 N · m · s and 172.12 N · m · s, respectively. So the desired control method consumes about 3% and 7% more than those from the discounted cost method for the acceleration and torque, respectively.

Fig. 9.

	Figure Option ViewDownloadNew Window
	Fig. 9. (color online) Control accelerations.

Fig. 10.

	Figure Option ViewDownloadNew Window
	Fig. 10. (color online) Control torques.

Figure 11(a) shows the trajectory of the chaser’s feature point with the desired control method in the x–y plane of the LVLH frame. The arrow points to the chaser’s feature point, and the direction is parallel to the feature point’s position vector in the body-fixed frame. Figure 11(b) shows the feature point’s trajectory in the x–z plane.

Fig. 11.

	Figure Option ViewDownloadNew Window
	Fig. 11. (color online) Trajectory of chaser’s feature point in planes of (a) x–y and (b) x–z.

In this paper, the problem of close-range approaching to a tumbling target is investigated. Firstly, a kinematically coupled 6-DOF model describing the relative motion of the spacecraft’s feature component is set up in the reference LVLH frame. The model can be used for the tumbling target and the maneuverable chaser individually. Then the proximity of the chaser for the tumbling target is defined as an optimal tracking problem. By adding a desired control signal into the index function, a feedback-feedforward tracking controller is designed. It is shown that the proposed controller is the derivative of the SDRE regulator, and has the suboptimal property. Compared with the existing discounted cost method, the proposed controller has a little computation advantage since the feedforward term is computed from the algebraic equation, not the Lyapunov equation. Numerical simulations show that the two methods can both perform the position and attitude tracking simultaneously for the tumbling target, and with a little more fuel consumption(less than 10%), the proposed method achieves evidently higher tracking accuracy than the discounted cost method.